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The numerical evaluation of an individual Bessel or Hankel function of large order and large 
argument is a notoriously problematic issue in physics. Recurrence relations are inefficient when 
an individual function of high order and argument is to be evaluated. The coefficients in the 
well-known uniform asymptotic expansions have a complex mathematical structure which involves 
Airy functions. For Bessel and Hankel functions, we present an adapted algorithm which relies on a 
combination of three methods: (i) numerical evaluation of Debye polynomials, (ii) calculation of Airy 
functions with special emphasis on their Stokes lines, and (iii) resummation of the entire uniform 
asymptotic expansion of the Bessel and Hankel functions by nonlinear sequence transformations. 

In general, for an evaluation of a special function, we advocate the use of nonlinear sequence 
transformations in order to bridge the gap between the asymptotic expansion for large argument 
and the Taylor expansion for small argument ("principle of asymptotic overlap"). This general 
principle needs to be strongly adapted to the current case, taking into account the complex phase of 
the argument. Combining the indicated techniques, we observe that it possible to extend the range 
of applicability of existing algorithms. Numerical examples and reference values are given. 

PACS numbers: 02.60.-x, 44.05. +c, 02.70.-c, 31.15.-p 



I. INTRODUCTION 



Bessel, Hankel and Airy functions constitute some of 
the most important special functions used in theoretical 
physics, and their calculation is notoriously problematic 
for extreme ranges of argument and order, even if their 
mathematical definition is straightforward. Especially, it 
should be noted that recurrence relations are useful for 
arrays of Bessel or Hankel functions when, for given ar- 
gument, all functions up to a maximum order are needed. 
However, recurrence relations cannot be used to good ef- 
fect if an individual function of large argument and order 
needs to be evaluated. 

Notably, Bessel, Hankel and Airy functions occur in 
the multipole decompositions of various operators in elec- 
trodynamics; these arc known to be slowly convergent 
decompositions in many cases. A lot of work has been 
invested into the development of asymptotic expansions 
which may be used in the calculation of the special func- 
tions. After the famous paper of Debye which used 
a saddle point expansion, a first treatise of the notori- 
ously problematic case of equal order and argument ap- 
peared in Ref. 0- Further historical papers, where the 
theory was refined, can be found in Refs. [3l-[lG|. The 
standard textbook [ll| contains a collection of very use- 
ful formulas. The known asymptotic formulas for large 
order (at fixed argument) are given in reference vol- 
umes (e.g., Refs. |12| - [l4| '). and they can be used, to- 
gether with asymptotic formulas for other Green func- 
tions [l^ , for the calculation of the properties of bound 
electrons. Indeed, Bessel, Hankel and Airy functions be- 
long to the most important special functions used in the- 



oretical physics; the much revived interest in these is also 
manifest in a recent monograph on Airy functions [l6j . 
The renewed interest in the theory of special functions 
is also manifest in a number of other recent books and 
review articles [l7H2(3l. 



In a marvelous tour de force, Olver |2ll . |22| | has derived 
uniform asymptotic expansions which hold for large order 
of Bessel and Hankel functions, uniformly in the complex 
plane of the argument variable (for arguments z with 
a complex phase |arg(z)| < tt — e). The derivation is 
based on the general asymptotic properties of solutions 
of second-order differential equations. The findings are 
summarized in Chapter 10 of the textbook [23| which is 
usually more accessible than the orig inal references (see 
also Chapter 8 of the recent Ref. |18j). 

A central question surrounding the use of the uniform 
asymptotic expansions has been their practical applica- 
bility to the calculation of Bessel and Hankel functions. 
This question is important because the expansions, while 
uniformly applicable in the complex plane, have a compli- 
cated mathematical structure, and because they involve 
Airy functions whose numerical evaluation is eventually 
required for arbitrary magnitude and complex phase of 
the argument. Despite considerable and perhaps justified 
doubts regarding their usefulness for numerical calcula- 
tions, the uniform asymptotic expansions [2ll[22j seem to 
be the most powerful ones available for the calculation of 
an individual Bessel or Hankel function of large argument 
and order. The aim of the current article is to show that 
their domain of usefulness can be drastically enhanced if 
they are combined with the "principle of asymptotic over- 
lap" that makes it possible to join asymptotic regions for 
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large argument with regions of small argument via the 
use of a nonlinear sequence transformation in overlap- 
ping regions [see also Section 2.4.1 of Ref. p3|]. 

The importance of the numerical evaluation of an indi- 
vidual Bcsscl functions for physics is highlighted by the 
substantial work devoted to the development of asymp- 
totic expansions and numerical algorithms |25l - [35| . If 
one aims to develop the algorithms based on the uni- 
form asymptotic expansions [2l], H^] , one first needs to 
evaluate Airy Ai and Bi functions. For large modulus 
of the argument and variable complex phase, their be- 
havior is characterized by a Stokes phenomenon. Along 
the Stokes lines, i.e., along specific values of the complex 
phase of the argument of the Airy functions, the relative 
magnitude the contribution of the different saddle points 
changes. A numerical algorithm for the Airy functions 
has been described in [36[. It relies on a separation into 
the region of large argument, where an asymptotic expan- 
sion is applied, and the region of small argument, where 
a differential equation is being integrated. Here, we ad- 
vocate the use of a nonlinear sequence transformation in 
order to bridge the gap between the asymptotic regime of 
large argument, and the regime of small argument where 
a power series can be used. The asymptotic expansion 
used for large argument has to be be adapted according 
to the complex phase of the argument. As a byproduct of 
our analysis, we derive some higher terms in the analytic 
expansion at the exact turning point v — z. 

The paper is organized as follows: We first recall ba- 
sic formulas in Section |lll Numerical calculations are 
described in Section Hill Analytic properties at the turn- 
ing point V = z are calculated in Section IIVI Finally, 
conclusions are drawn in Section [V] The Appendix |X] is 
devoted to a general discussion about the saddle points 
in the complex plane related to the Bessel functions, and 
about the possibility of constructing an alternative algo- 
rithm. 



II. BASIC FORMULAS 

For complex argument z with Re(z) > 0, the evalua- 
tion of Bessel J functions can be traced to the evaluation 
of integrals of the form [see Eq. (10.9.6) of Ref. [l3| ] 

1 r 

J^{z) = - cos [z sm{e) -iy0]de 
Jo 

_sm{un)_ /■°°exp[-z sinh(0) -zy0]d0, (1) 
7^ Jo 

where the order f of the Bessel function is not necessarily 
an integer and | arg(z)| < 7r/2. Of particular interest are 
the Bessel J functions, as they are regular at the origin 
for positive integer ly. For integer v, the second term in 
the definition of J according to Eq. ([T]) vanishes. All of 



the definitions used here for Bessel functions, and spher- 
ical Bessel functions, are contained in Chaps. 9 and 10 
of Ref. [T^]. Indeed, these and many of the asymptotic 
formulas used in the following are also included in the 
modernized handbook [l^ [ij]. Reference [l^l is now 
somewhat outdated but still the standard classic refer- 
ence on the matter. 

For half-integer ly, the Bessel J functions are related 
to spherical Bessel functions according to the formula {i 
is an integer) 



(2) 



This relation is given in Eq. (10.47.3) of Ref. For 
Rc(z) > 0, the Bessel Y function is defined as [see 
Eq. (10.9.7) of Ref. [H] 

y^(z) = - / sin [z sin{e) - u 6] AO 
I" Jo 

{e'^* +e-'^*cos(z/7r)} e"^ '^'"''(*M6' , (3) 



for I arg(z)| < 7r/2. The spherical y function is defined as 



Yi+i/2{z) . 



(4) 



This relation is given in Eq. (10.47.4) of Ref. [l3| . 
The Hankel functions are defined as [see Eqs. (9.1.3) 
and (9.1.4) of Ref. [13] 



hI^^z) = Mz) + iY,{z). 

hI^\z) ^ J„iz) - iV^iz) . 
h['\z)=Mz)+iyeiz), 
^''e^\z) = it{z) -iye{z) . 



(5a) 
(5b) 
(5c) 
(5d) 



The definitions of the spherical Hankel functions are 
given in Chap. 10.1.1 of Ref. 

The integral representations Q and (jS]) are valid for 
Rc(z) > 0. For purely imaginary z, we may use a def- 
inition in terms of the modified Bessel functions I,^{x) 
and K^{x), see Eqs. (|A4p - (jA7p . Below, we describe a 
numerical algorithm with the notion < arg(z) < tt in 
mind. Arguments z with — tt < arg(z) < are treated 
by the transformation 2 — > z' = z exp(i tt) (so that the 
transformation z — )■ z' does not leave the first Riemann 
sheet). The conversion formulas can be derived based on 
Eqs. (9.1.35)— (9.1.39) of Ref. [H and read 

J^{z)=e-''''' J^ize'''), -TT < arg(z) < , (6a) 

Y^{z) ^ Y^{z c'"") ~ 2i cos{iy tt) J^izc'"") , (6b) 

i/W(^) -c''^"i/W(ze'")-f 2e-'''" J,(ze'"), (6c) 

^ _ gi . vr e' ") . (6d) 



3 



(Many of the asymptotic formulas given here are also 
included in the modernized handbook [l^, but for the 
time being, wc prefer to refer to equation references in 
the somewhat outdated, but standard classic Ref. [l^.) 
Alternatively, one may use direct complex conjugation, 

Mz) = iMz*))* , Y,{z) ^ {Y,{z*))* , (7a) 

i7W(z)= (i/(2)(z*))*,ij(^)(z)=(i7W(z*))*, 

(7b) 

where z* is the complex conjugate of z. Using a combi- 
nation of the formulas © and (O, we could in principle 
restrict the range of complex phases of the arguments 
to the first quadrant of the complex z plane. However, 
as evident from Eqs. (US]), (UHl), ^ and below, wc 
would still need to evaluate the Airy Ai and Bi functions 
in the entire complex plane, even if we restrict z to the 
first quadrant in the complex z plane (and the former 
constitutes the main computational challenge) . A simple 
restriction to the upper half of the complex z plane thus 
seems to be most effective. 

Without loss of generality, we restrict our attention to 
the case > in the following. For < 0, the conversion 
formulas are as follows, 

H^'hz) = e'^'' Hl^\z) , (8a) 

H^^liz) ^ c-'^'^ Hl^\z) , (8b) 

J-u{z) = cos(7ri/) Ji,{z) — sin(7rt^) Y'^(z) , (8c) 

y_^(z) = sin(7ri^) J^(z) + cos(7ri/) y^(z) . (8d) 

The conversion matrix for the Bcssel functions J and 
Y has the same structure as a rotation matrix for an 
angle ttv. For the Hankel functions, the above formu- 
las ([Sal) and (l8b| can be found in Eqs. (9.1.5) and (9.1.6) 
of Ref. [13. 

In principle, one might speculate that the above inte- 
gral representations ([T]) and ([3]) should be sufficient in 
order to numerically evaluate an individual Bessel func- 
tion. However, the numerical difhculties for large u arc 
nearly insurmountable in view of apparent numerical os- 
cillations of the integrand. While one can investigate 
complex integration contours with the notion of adopt- 
ing a steepest descent method (see Appendix VK\ . these 
representations do not immediately lead to a uniformly 
applicable algorithm, either. 

Finally, let us recall the basic asymptotic properties of 
Bessel J and Y functions for > 0. Only the Bessel J 
functions is regular at the origin, and we have 

y,(z)^ -^0) , 1^1 ^0. (9b) 



The literature on Bessel functions is manifold. A very 
useful reference is the standard treatise [ll|. In Chap- 
ter 10 of Ref. [i^, basic asymptotic expansions and prop- 
erties of Bessel J and Y functions, and of their deriva- 
tives, are reviewed and explained very clearly. 

One is often faced with the problem of calculating 
strings of Bessel functions whose indices differ by inte- 
gers [37H39j . Recursive algorithms based on the relations 

Ju-i{x) + Ju+i{x) = — Jv{x) , (10a) 

X 

r,_i(a;) + r,+i(x) = — Y,{x) , (10b) 

X 

can be very effective, as explained in Section 10.5 on 
p. 452 of Ref. [l^. For Bessel J functions, one starts a 
three-term downward recursion in v with two essentially 
arbitrary starting values for Ji^j^iix) and Jv{x) at high 
V and continues to calculate Ji,-i{x) until the order of 
the Bcssel function becomes zero. One can then either 
calculate Jo {x) explicitly and use the recurrence relation 
upwards (filling the array of Bessel functions), or fix the 
normalization of all calculated Bessel functions by a nor- 
malization condition |40| - |43 | (see also Chap. 10.10.5 of 
Ref. [13 )■ In Ref. (42| . the computational aspects of 
three-term recursion relations have been discussed with 
a special emphasis on their numerical stability. Here, we 
are dealing with a different problem, namely, the evalu- 
ation of an individual Bessel function of high order and 
argument Ju{x) and Yu{x), without recourse to any re- 
currence relation in v. 



III. SUMMATION OF THE UNIFORM 
ASYMPTOTICS 

A. Uniform asymptotic expansions 

The task in the current investigation is to calculate the 
functions 

J^{z^vy), Jl{z = vy), (11a) 

Y,{z = yv), Yliz^yy) (lib) 

for complex argument — tt < arg(z) < tt, and real v > 0. 
In a numerical code, it is sufficient to treat the complex 
phase range < arg(z) < tt. The range — tt < arg(z) < 
is covered by Eqs. © and ([7]). Because we arc us- 
ing asymptotic expansions valid for large v, we also as- 
sume that V > 50. For i/ < 50, one may use Miller's 
method [iol - IZ^ . A brief digression on this algorithm can 
also be found in Chap. 10.10.5 of Ref. [l2|. The case 
1/ < —50 then is covered by Eq. ([8]). The parameteriza- 
tion z ^ vy is useful to identify the notoriously prob- 
lematic region near y ~ 1. It is also being used below in 
Appendix [X] 
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We first have to recall the uniform asymptotics of 
the Bessel functions (see Refs. [21], [22j), These are also 
listed in Eqs. (9.3.35), (9.3.36), (9.3.43V and (9.3.44) of 
Ref. [HI, and in Chapter 10 of Ref. [2|. A brief red- 
erivation is given in Ref. [29| . The uniform asymptotics 
for the Bessel J function are given by 



4C 



1/4 



Ai(^^/^C) MO 



1/3 



,2k 



k=0 



1/5/3 



jy2k 

k=0 



(12) 



This asymptotic formula is valid for v — > ±00 and 
^^'&{y) — where e is an arbitrarily small positive 
number. We denote the Airy function of the first kind 
as Ai. 



For y > 0, the ( variable is defined as 



2/3 



2/3 



0<2/<l, 

2/3 



— 1 — arccos I - 

y 



2/3 



< 0. 



y > 1 



>o, 

(13a) 
(13b) 



The calculation of C for complex y relies on the formula 



|<./^.,„ i±v^ -/r^ („, 



where the branches take their principal values when z € 
(0,1) and C € (0,oo) and ( is continuous elsewhere, as 
described in Chapter 10.1 of Ref. [23| . 



The Ofc and bk coefficients entering Eqs. (IT2]) — ((20 
read 



2k 



ak 



(15a) 



s=0 



2fe+l 



(15b) 



s=0 



formulas read 



3! 144^ 



m^2s+l 
m odd 



r(3s + ^) 
9V^r(2sTi) 



1 



s!l44 
6s + 1 



(2s + 1) (2s + 3) ••• (6s-l) 
2r(3s+ I) 



(16a) 



6s- 1 



A. = - 



9''0F(6s- l)r(2s+ 1) 



-. (16b) 



The Debye polynomials fulfill uo{t) = 1 and arc otherwise 
defined recursively as 

Uk+i{t) = i (1 - t') 4 W + 1 t dt' (1 - 5t'2) Uk{t') . 

o Jo 

(17) 

For polynomials, the operations of differentiation and in- 
tegration can be represented by simple multiplication op- 
erations acting on a coefficient matrix. This is due to the 
trivial identity dx"'/dx — nx"~^, applied to integer n. 
On a computer system, it is thus possible to evaluate the 
coefficients of, say, the polynomial coefficients for the first 
few hundred Debye polynomials and to use them in order 
to evaluate the ak{C) and &/c(C) coefficients for given (. 

The asymptotic expansion (jl2p obviously is an expan- 
sion for large and it is valid even in the problematic 
region y « 1. We are now in the position to give the 
corresponding formula for the Y function, which involves 
the Airy function of the second kind Bi and its derivative. 



^ < jy 

k=0 



2 k 



,2k 



(18) 



k=0 



The uniform asymptotic expansion of the Hankel H^^^ 
function is given by 



i7W(..)^2e"-/3 (^) 

f Ai(e2W3 j,2/3 ^) - 



1/4 



1/1/3 ^ J. 

fe=0 



2k 



,2.i/3 Ai'(e2-'/3 ^2/3^ y MO 1 ,^g. 



„5/3 



k=0 



They involve the Debye u polynomials, and coefficients 
/is and Xg which need to be defined. The corresponding 



According to Eq. (|5al) , the uniform asymptotic expansion 
for ij'^-* is obtained by changing the sign of the imaginary 




Re(x) 




5 fsW 



(a) 





9b(x) 



FIG. 1: (Color online.) Figure (a) shows a contour plot of 
fA{x) = \Ai{x)\^^^ as a function of Re(a;) and Im(a;). The 
exponent i is introduced in order to prevent "overflow" of 
the plotted function near the boundaries of the considered 
range of arguments. Figure (b) shows a contour plot of 
I exp( — I a;^''^)!^^^ as a function of Re(a::) and Im(a;). For 
X < Q, i.e. on the negative real axis, the modulus is unity 
(Stokes line) . The zeros of Ai(a:) give rise to the visible "bump 
holes" on the negative real axis in panel (a). Except for the re- 
gion near arg(a;) = ±7r, the Airy Ai function can be described 
using a single, uniform asymptotic formula A{x) as defined in 
Eq. (|24a|l . which is proportional to | exp(-| x^^'^)\^/^ . 



FIG. 2: (Color online.) Figure (a) shows aplot of |Bi(a;)|^'''' as 
a function of Re(a;) and Im(2;). From the contour plot, it is ev- 
ident that a simple exponential of the form | exp(±| x'^^'^)\^^^ 
cannot possibly describe the asymptotic behavior of the 
Airy Bi function. The zeros of Bi(a;) give rise to visible 
"bump holes" on the negative real axis and along the lines 
arg(a::) = ±7r/3 in panel (a). Figure (b) shows a contour plot 
of |/(a:)| where f{x) = \exp(^ x^/^)\^/'^ for |arg(s:)| < 7r/3 
and f{x) = |exp(-|a;3''2)|i''® for tt/3 < |arg(x)| < tv, rep- 



resenting the Stokes line behavior near arg(a;) 
arg(a;) = tt [see Eq. (f27ll ]. 



±7r/3 and 



unit, 



B. Evaluation of the Airy function 



4C 



1/4 



Ai(e-2-/3 ^2/3 - akiO 



A/3 



k=Q 



,2k 



g-2W3 Ai'(e-2-/3^V3^ g5fc(C) 



.5/3 



,2k 



(20) 



k=0 



The corresponding equations for the derivatives of 
the Besscl and Hankel functions can be found in 
Eqs. (9.3.43), (9.3.44) and (9.3.45) of Ref. [13. Oth- 
erwise, the derivatives are also accessible via the formula 

^^Jl{z) ^ Jl[z) = \ {J,^,{z) - J,+,{z)) , (21) 

where J stands for J, Y, i/^^' or H^'^\ 



We now discuss the principle of asymptotic overlap in 
the evaluation of the Airy functions of the first and sec- 
ond kind. In contrast to the usual notation, we denote 
the complex argument of the Airy functions as x, in or- 
der to distinguish it from the argument z of the Bessel 
and Hankel functions. For x — )■ 0, we use the expansion. 



Ai(a;) 



V— ^ 

f-^ k\ T (k 

k=0 V 



„3 k 



OO 



,3 fc-l-1 



(22a) 



For the derivative of the Airy function, we have 
1 



Ai'(x)= -J2 



2-2fc-3 



„3 k 



^^^fc!r(fc + f) 



^3fc+2 



(22b) 
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The Airy Bi function is given by 



functions at large argument, 



Bi(a;) 



oo 2k-'i 



^3fc+l 



(22c) 



and its derivative by 



Bi'(a;) = £ 



oo 

„3fe , 



3-2fc-g 



^3 fc+2 



-^fc!r(fc + i) ^^^fc!r(fc + f) 



(22d) 



The convergence radius of the expansions ([22|) is actu- 
ally infinite, because of the Gamma functions in the de- 
nominator. Still, they are faced with numerical prob- 
lems for large negative x. In order to illustrate this 
fact, we draw an analogy to the (likewise convergent) 
expansion exp(— x) = I ^ which also in- 

volves a power in the numerator and a factorial in the 
denominator. Being absolutely convergent, this expan- 
sion is numerically useless for the evaluation of exp(— a;) 
at X = 5000, because the largest term in the series is 
of order exp(5000), whereas the entire sum amounts to 
exp(— 5000) ~ 3.37 x 10~^^^^. If the series were summed 
term by term, one would incur a numerical loss of more 
than 4000 decimals. Other methods (asymptotic expan- 
sions for large argument) therefore have to be pursued. 
Figures [T] and O show the Stokes phenomenon. From 
Fig. [21 we infer that the Airy Bi integral is exponen- 
tially growing for |a;| — >■ oo, along the complex directions 
arg(a;) = and arg(a;) = ±27r/3, with Stokes lines at 
arg(a;) = ±7r/3 and arg(a:;) = tt. Therefore, the Bi in- 
tegral cannot be represented by a simple asymptotic di- 
vergent series but must be the sum of two. This can be 
justified on the basis of saddle point considerations [2]| . 



The expansions ([22)1 are valid for a; 0. Comple- 
menting these expansions, for x — >■ oo, one is interested 
in suitable asymptotic expansions of the Airy functions. 
To this end, it is useful to relate the Airy Ai and Bi func- 
tions to a modified Bessel functions K and / of order ±-i. 



Ai(a;) 
Bi(a;) 



,3/2 



Mir 



3/2 



(23a) 
(23b) 



Based on the relationship of the Airy functions to modi- 
fied Bessel functions, we infer that the following asymp- 
totic scries are relevant for the calculation of the Bessel 



exp (-| x^/^) 
7r3/2 .t1/4 



- (-3)'-T(fc + i)r(fc + t) 

22k+2 ^3fc/2 ' J 



B{x) 



k=0 

exp (I x^/^) 
7r3/2xi/4 



22k+2 ^[ 2.3fe/2 



(24b) 



C{x) 



k=0 

a;i/4 exp(-|x3/2) 

,v (Vr(fc-i)r(fc-f-i) 

^ 22k+2 ^3k/2 ' 



fe=0 



whereas D{x) is the only function that carries an overall 
minus sign in the pref actor, 



Dix) 



xi/4 exp (f x3/2) 



7,3/2 

^ 3"r(fc^^) r(fc + |) 

^ 22fc-|-2 r^3k/2 • 



fe=0 



These series are generalized hypergeometric series of the 
2^0 form that diverge for every nonzero argument .x^3/2 
unless they terminate. Their asymptotic relation to the 
Airy functions is clarified below. Indeed, the asymptotic 
expansions of Ai and Bi and of their derivatives can be 
related to the series A{x), B{x), C(x) and D{x). For 
|x| — >■ oo, we have the divergent asymptotic expansion as 
a function of the complex phase, 

Ai(x) ^ A{x) , 

|x| — )- oo , I arg(x)| < TT — e , (25a) 

Ai(x) ^ A{x) +iB{x) , 

|x| — oo , n — e < arg(x) < tt , (25b) 

Ai(x) - A{x) -iB{x) , 

|x| — J- oo , — TT < arg(x) < — TT -I- e , (25c) 

For the derivative of the Ai function, we have 

Ai'(x) - C(x) , 

|x| — )■ oo , I arg(x)| < TT — e , (26a) 

Ai'(x) - C(x) + ii:i(x) 

|x| — oo , TT — e < arg(x) < tt , (26b) 

Ai'(x) - C{x)~iDix) 

|x| — > oo , — TT < arg(x) < — tt 4- e , (26c) 

for arbitrarily small e. The role of the e parameter in 
these expressions can be clarified as follows. From the 
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formulas, it is evident that arg(a;) — >■ ±7r (from below or 
above, respectively), there is an admixture of the asymp- 
totic B and D series to the A and C series. The magni- 
tude of that admixture is related to the magnitude of 
For large |a;|, we can choose e to be small. For finite e, one 
has to compare the magnitude of the two contributions. 
If the subdominant saddle point is numerically signifi- 
cant, one has to add its contributions (see also Fig. [3] 
below). 

The dominant contributions to the asymptotic expan- 
sions for the Airy Bi functions are given by 



Bi(a;) - 2B{x), 



^ oo , I arg(a;)| < 3 ^ > 



(27a) 



Bi(a;) - 2B{x) +iA{x) , 



\x\^oo, ^ -e < arg(a;) < ^ -He, (27b) 



Bi(a;) r^iA{x), 



\x\ 00 , — + e < arg(a;) < tt — e , (27c) 



Bi(a;) - Bix)+iA{x) , 

|x| — > 00 , TT — e < arg(a;) < tt . 



(27d) 



For arg(x) < 0, the corresponding conditions are ob- 
tained by complex conjugation, 

Bi(a;) - 2B{x) -iA{x) , 

\x\~^oo, -| - e < arg(a;) < + e , (27e) 

Bi(a;) iA{x), 

\x\ 00 , — — + e < arg(a;) < — tt + e , 

Bi(a;) ~ B{x) - iA{x) , 

|a::|->oo, -tt < arg(x) < -tt + e . (27f) 

The asymptotics for the derivative of the Airy Bi function 
are obtained by making the following replacements in 
Eq. dsn): Bi ^ Bi', A^C, and B ^ D. 

The "principal of asymptotic overlap" which has been 
employed for a number of special functions in Ref. (2^ 
would now call for an application of the power se- 
ries ([22a|) . ((22bl) . ([22cl) and ([22dl) for smaU \x\ < Ro, 
and for the use of the asymptotic expansions (ps)) . ([25]) 
and ([271) for |a;| > with the region Rq < \x\ < Ri 
being bridged by the application of a nonlinear sequence 
transformation to the asymptotic expansions ([25|) . ([26| . 
and (|27p . However, the situation is more complicated in 
reality, and it is impossible to choose the parameters Rq 
and i?i uniformly in the complex plane, independent of 
the complex phase of the argument x. 

A possible scheme which is sufficient to obtain at least 
20 digits of accuracy for all x is outlined in Fig. [3] The 



designation "POWER" indicates the use of the power 
series (|22a|) . (|22bp . ([22cl) and (|22dp . as appropriate, the 
designation "ASYMP" indicates the application of the 
asymptotic expansions (PS)) . (pS)) . and ([TT]) . where the 
series A{x) and B{x) are given in Eqs. (|24a[) and (j24bp . 
Finally, the designation "TRAFO" needs to be explained: 
it denotes the application of a nonlinear sequence trans- 
formation to the asymptotic series, with the aim of ex- 
tending their validity to lower modulus of the argument 
than what they would otherwise be valid for. 

A rather thorough discussion of an application of se- 
quence transformation to a nontrivial problem has been 
given in Section 2.4.2 of Ref. 



24 



in the context of the 
relativistic Green function for the hydrogen atom. Es- 
sentially, sequence transformations are generalizations of 
Fade approximations that fulfill accuracy-through-order 
relations, i.e., they constitute rational functions that re- 
produce the first few terms of a given input series when 
expanded back in powers of the argument. By reformu- 
lating the sequence transformation in terms of optimized 
remainder estimates [i^ , one can enhance the rate of con- 
vergence as compared to Fade approximations. A more 
thorough discussion of sequence transformation would go 
beyond the scope of the current article. The sequence 
transformation used in the current article is defined as 
follows. Let Sn denote the nth partial sum of a series 



k=0 



ak 



1,2, 



(28) 



where the a„ are the terms in the series to be summed, 
and in our case, the terms in the asymptotic expan- 
sions ([24)1 . Let the difference operator A be explained 
as 



(29) 



As described in Refs. [24, l43|, |4J] (see also §3.9 of 
Ref. [13II). in many cases the following sequence trans- 
formation (Weniger transformation), 



J=0 



il3 + n + j)k- 



{l3 + n + k)k-i Asn+j 



k 

E 

j=0 



i-iy 



{l3 + n + j)k^ 



1 



{P + n + k)k-i As„+j 

. (30) 

leads to an accelerated convergence (or summation in 
the case of divergence) of the input sequence {sn}J5^o 
of partial sums of the series of the a^. Here, /3 is a 
shift parameter which we choose as /3 — I (as given 
in Refs. [13, 13 )■ The starting order for the transfor- 
mation in Eq. ([50)) is n which we choose to be rt = 0. 
Then, 5[°''(/3,so) denotes the fcth order Weniger trans- 
form of the input sequence {sn}^o- The central idea in 
the construction of the sequence transformation ([50]) is 
the expansion of the remainder term in an inverse facto- 
rial scries, as explained in Ref. [43j . The use of recursion 
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relations described in Section 8.3 of Ref. [43| is impera- 
tive in order to ensure the computational efhciency and 
numerical stability in the computations of the Weniger 
transforms. In Refs. [13,111,14^, it has been established 
that the Weniger transformation is very powerful at sum- 
ming the factorially divergent series that result from the 
asymptotic expansions of special functions in terms of 
hypergcomctric series. 

We find, by comparison to calculations with extended 
arithmetic using computer algebra systems [4^ and by 
comparison of different methods along the separating do- 
mains indicated in Fig. |3l that naive termination criteria 
are sufficient in order to reach a prescribed accuracy of 20 
decimals in the entire complex plane, for the Airy func- 
tions. By a naive termination criterion, we mean that 
the summation of terms in the power series, or the sum- 
mation of terms in the asymptotic series, is terminated 
when the next higher-order term divided by the most re- 
cently calculated partial sum is a factor 100 less than the 
prescribed accuracy for the evaluation of the Airy func- 
tion. Likewise, the calculation of subsequent higher-order 
transforms d'j^^ of the Weniger transform (pO)) is termi- 
nated when the apparent convergence of three consecu- 
tive transform (maximum absolute value of the relative 
difference of transforms A;, A; -1-1 and k + 2) is smaller than 
the prescribed accuracy by a factor of 100. 

The appropriate evaluation method for given complex 
argument x is indicated in Fig. |31 as a function of Re(a;) 
and Im(x). The derivatives Ai'(x) and Bi'(.T) are evalu- 
ated using the same schemes as Ai(x) and Bi(a;), respec- 
tively, but (for the asymptotic regime of large argument) 
with the replacements A — )- C, and B ^ D [see Eq. ([M)) ]. 
We should also clarify the exact formulas for some of the 
separating curves in Figs. |3Ja) and EJb). The curved 
separating line between "TRAFO" and "POWER" in 
Fig. IHa) follows the formula 



kl < 5- 



15 



arg(a;)|, |arg(x)|< 



2tt 



(31) 



and it meets the outer asymptotic region "ASYMP" at 
|arg(x)| = ^ and = 15. In Fig. [HJa), the transition 
from the asymptotic scries A to the asymptotic series 
A + iB takes place at arg(x) = In Fig.[3Ib), for the 

Bi integral, the power series is used for \x\ < 5 uniformly 
for any complex phase of cc. For |x| < 15, the power series 
is also used, but only for | arg(a;)| < 7r/3. The transition 
from the asymptotic series 2 B to the asymptotic series 
2 B + lA takes place for \x\ > 15, and arg(a;) = ±7r/6. 



C. Summation of uniform asymptotics 



The uniform asymptotic formulas ([11]), ([TH]), 
and (|20p arc expansions for large v for argument z — uy 
of the Bcssel functions. They can be written in a form 





FIG. 3: (Color online.) Algorithms used for the Airy Ai func- 
tion [panel (a)] and for the Airy Bi function [panel(b)]. The 
explanation is in the text. 



that is amenable to the application of the convergence ac- 
celeration transformation (|30p . as follows. We consider 
as an example Eq. (|12p . Then, 



k=Q 



k=0 



(32a) 



ak 



4C \'/' (Ai{,.y^O a,{C) 



A/3 



,2k 



y5/3 



,2k 



(32b) 



Here, the Airy functions need to be evaluated only once 
as they do not depend on k. Treating the Airy function 
term and the derivative term together (i.e., as a single 
term dk), we empirically observe better convergence in 
many cases. Written in the form (|32|) . the convergence 
of the uniform asymptotic for the Bessel J function can 
be accelerated by calculating the transforms (5["^(/3, s„) 
where again, we choose n = in the current investiga- 
tion. We empirically observe that the use of the con- 
vergence accelerator does not induce any stability issue 
even in cases {v ^ z) where the apparent convergence of 
the asymptotic series would be sufficient to calculate the 
Bessel J function to the required accuracy. 

A special case still necessitates a modification of the 
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algorithm. Namely, for large and almost equal v and 
z, the confluence of the two saddle points defining the 
Bessel function (sec Appendix lAl below) leads to numer- 
ically prohibitive cancellations in evaluating the series of 
the afc(C) and 6fc(C) coefficients. It is impossible to over- 
come these difficulties with the necessarily finite precision 
of computer systems in the limit z — > i/; indeed, as de- 
scribed below in Section IIVI the asymptotic expansions 
for large v at exact equality v = z are obtained after 
the cancellation of a number of divergent terms in e after 
setting V ^ z + e. 

Eventually, we find that the region in very close vicin- 
ity of the turning point v = z can only be overcome by 
an explicit use of the recurrence relation 



For an example with arg(z) = 7r/3, we obtain 



(33) 



with the notion of displacing v from z, where again 
stands for J, Y, i?^^^ or H^^\ We find that numerical 
cancellations are most severe when v is slightly lower than 
Re(2;). So, for Re(z) — 21 < ly < Re(z) -I- 2, we therefore 
express J71/(z) as a function of j7i/+24(2) and Ji,+25{z) by 
repeated use of the recurrence relation ([55)) . Indeed, by 
repeated application of Eq. p3p one can express Jy{z) 
as 



Ju{z) = f2i{,V-, z) Ju+2i{z) + f2b{v, z) Jy+2b{z) 



(34) 



where /24 and /25 are somewhat lengthy rational func- 
tions of their arguments. Their explicit form can easily 
be obtained using computer algebra [46| . 



The recurrence relation p4|) is applied in the direction 
of increasing order v of the Bessel function, where it may 
be be unstable [indeed, Ju{z) — >■ and Y^{z) — >• cx) for 
1/ — >■ oo at constant z\. One might thus expect that the 
shift (p4| could induce numerical instability. However, 
because the shift is applied in the region v z, where 
the Bessel function is not yet in its asymptotic regime for 
large 3> z, the numerical cancellations are only minor 
and do not constitute a matter of concern. 



D. Numerical examples 



In order to demonstrate the power of the algorithms 
proposed in the current article, we now turn our attention 
to a number of numerical example cases. We first present 
two evaluations for non-integer order and argument of 
Bessel functions, in the range ;/ w z, which result in 

^5000000 2(5 000 000.1) = 2.614 463 954 691 926 x 10"^, 

(35) 

and in 



5 000 000 



.2(5000000.1) 



-4.533 251 771 400 041 x 10"^ . 

(36) 



H, 



5000000.2(5000000.1 exp(iV3)) 

= -6.120 398 939 598 734 x io-954990 

- i 1.992 559 471 616 042 x 10-954989 ^ 



(37) 



While the result is nearly zero, numerical evaluations of 
this kind are needed because other terms in angular mo- 
mentum expansions in quantum electrodynamics lead to 
extremely slowly convergent series (47l . l48j due to an in- 
terplay of increasing and decreasing terms (as the angular 
momentum quantum number is increased). Another ex- 
ample of explicit evaluation of the Hankcl function for 
extremely large order and argument in the turning point 
regime 1/ « z reads. 



(1) 



6000 000.2(6 000 000.7) 

= 2.467 848 322 382 092 x 10"^ 
- i 4.252 887 224 934 845 x 10" 



(38) 



and 



(39) 



6 000 000.2(6 000 000.7) = 

2.467 848 322 382 092 X 10"^ 

+ i 4.252 887 224 934 845 x 10"^ . 



For r « 1, it is instructive to verify the Bessel functions 
using the following sum rule, 

for r G (0, 1) and y > 0. The sum over £ constitutes 
a slowly convergent series whose can otherwise by accel- 
erated using the combined nonlinear-condensation trans- 
formation [49|. The sum rule directly follows from the 
angular momentum expansion of the Green function of 
the Helmholtz equation for the case of collinear argu- 
ments, as given in Chap. 9 of Ref. [s^l- It can also be 
derived as a reformulation of Eq. (10.60.3) of Ref. [isj . 
We have checked our values of the Bessel functions on 
the basis of this sum rule and plotted Bessel functions of 
high order in the turning point region (see Figs. H] and [S]). 



IV. ASYMPTOTICS AT THE TURNING POINT 



The direct application of the uniform asymptotics be- 
comes problematic when the argument and the order of 
a Bessel function are almost equal, because of numerical 
cancellations involved in evaluating the individual coeffi- 
cients in this case. We remember that for the case j/ w z, 
special methods have to be employed for the summation 
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0.002 
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19999000 


-0.001 
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FIG. 4: Plot of the Bessel function Ju{x) with v 
20, 000, 000.2 in the turning point region. 



0.002 




19999000 














-0.015 





FIG. 5: Plot of the Bessel function Y^{x) 
20, 000, 000.2 in the turning point region. 



with 



of the uniform asymptotic expansion (sec Section [ill Cp . 
This is because of highly significant numerical cancella- 
tions in the evaluation of the sums defining au and bk 
given in Eqs. pSap and (jl5bp . These numerical cancel- 
lations refiect the confluence of the two saddle discussed 
below in Appendix 1X1 and are present even if the sums in 
Eqs. ([T5a)) and ((T5b| are finite. 



Nevertheless it is possible to investigate the limit z/ — >■ z 
by an analytic expansion setting v = z -\- e. In the calcu- 
lation, the result is obtained only after the cancellation 
of a few divergent terms which are inverse powers of e. 
For exact equality i/ = z, it is possible to derive analytic 
approximations. These serve as an important test device 
for the numerical algorithms. The results given in for- 
mulas (9.3.31)"(9.3.34) of Ref. [l2| for lower order terms 
in the asymptotic expansions do not reach a sufficiently 
high order in order to be useful for the current investiga- 
tion, and we have thus calculated a few more terms in the 
asymptotic expansions of Jy{v) and Yu{v) for v ~^ oo. 



The general structure of the asymptotic expansions is 



as follows, 

J» = 

y» = 



jyl/3 2^ y2k ,,5/3 2^ T,2k ' 
31/2 



jy5/3 ^ , 

k=0 



00 

■a aj^ 

2 k 



j,l/3 ^ y 

k=0 



^5 / 3 jy2 k 



Jk_ _ a Ok 

,,2/3 ,ak ,,4/3 2^ ,,2k ' 



l/2/j — < y. 

k=a 



yi/3 L < y^ 

k=0 



31/25 



y2/3 L < y 

k=0 



E 



7fc 

2 k 



31/2 g Sk 

^4/3 Z-^ y2 k 
k=0 



The constants a and b are given by [IS 

21/3 



and 



32/3r(2/3) ' 
22/3 



(41a) 
(41b) 
(41c) 
(41d) 

(42) 
(43) 



3i/3r(i/3) 

Numerical values of the coefficients ak, 7fc, 5k were 
given in (unnumbered) equations following Eq. (9.3.34) 
of Rcf. [l2| . but only in numerical form. Using the for- 
malism outlined in Refs. [HI [13: Uli ; and computer alge- 
bra [4^ , it is relatively easy to evaluate these coefficients 
analytically. The first terms read 

1 151439 



as 



04 



1 CKi = Ck9 — 

225 ' 218295000 
887278009 
^ 2504935125000 ' 

1374085664813273149 
3633280647121125000000 ' 

1065024810026227256263721 



as = . (44) 

1540965154460247140625000000 ^ ' 

for the afc coefficients. For Uj with j = 6, 7, 8, the inte- 
gers in the numerator and denominator arc too long to 
be displayed, practically. The results read, up to 30 dec- 
imals, 

as = 0.00192 82196424877570138042 30112, 
ay = - 0.00762 60912 66562 73551 11507 07644 , 
as = 0.04059 16252 02439 02610 46110 96353 . (45) 
We obtain for the Pk: 

Pi 

1UZO(OU 

9597171184603 



/3o 
P2 



1 „ 1213 

70' '"^ 1023750 
16542537833 



,/33 



37743205500000 ' 25476663712500000 
53299328587804322691259 



/35 = 



91182706744837207500000000 ' 

70563256104582737796094772987 
~ 49341242187300033908437500000000' 



(46) 
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For Pj with j — 6, 7, 8, the results read, up to 30 decimals, 
/36 = 0.00506 84595 77410 25775 49192 90289 , 
/37 = - 0.02455 31387 44039 61140 77960 53629 , 
/38 = 0.15584 22313 83604 27406 01873 20654 . (47) 

The results for the 7^ read, 

23 604523 



70 = 1 , 71 



3150 ' 
2264850139339 



72 



' 644962500 ' 



73 = 

74 = 

75 = 



5095332742500000 ' 

160913976870912403 
^ 353106559055250000000 ' 

216363828773939104866579281 
266709417228648831937500000000 ' 



(48) 



For 7e, 77, and 73, we again give numerical results, up to 
30 decimals, 

76 = - 0.00222 15829 52781 90513 99808 25549 , 

77 = 0.00866 43390 74500 57903 06978 05128 , (49) 

78 = - 0.04561 48821 43058 95740 37593 20268 . 

Finally, we obtain for the ^fc. 



s - ^ A - 947 
'^''"5' '^^""346500' 

11192989 , 

02 = , 03 



(50) 



100443412440047 



18555075000' " 262141460831250000 
11007971229145235539 



(55 = - 



22905464949241875000000 

180026595127347603856424798473 



180354561678027325338750000000000 
The first 30 decimals of the results for ^6,7,8 read 

^6 = 0.00309 74954 54537 99792 65245 46715 , 

h = - 0.01341 95416 64085 43255 61825 95640 , 

(58 = 0.07735 88016 72766 55151 01508 73624 . (51) 

We have used these analytic formulas in our verification 
of numerical results obtained using the method described 
in Section Hill 

V. CONCLUSIONS 



We have described an algorithm for the evaluation of 
an individual Bessel or Hankel function for large order 
I/, and arbitrary complex argument z. The method relies 
on the use of the uniform asymptotic expansions given in 



Eqs. (IT2]), (HH), (HH) and These asymptotic expan- 

sions involve Ai and Bi functions, and their derivatives. 
Consequently, in passing and also as a prerequisite for our 
algorithm, we need to develop numerical code for the Airy 
functions valid in the entire complex plane. This is de- 
scribed in Section IlIIBi and the appropriate algorithms 
for the different regions in the complex plane are given 
in Fig. 131 Essentially, these represent an adapted imple- 
mentation of the "principle of asymptotic overlap" where 
a power series is used for small argument, an asymptotic 
expansion is used for large argument, and the region in 
between is bridged by a nonlinear sequence transforma- 
tion [Weniger transformation, see Eq. (|30l) ]. 



As evident from Fig. [31 this principle needs to be 
adapted for the Airy functions, and the regions separat- 
ing the algorithms depend on the complex phase of the 
argument. Using the Airy function algorithm and the 
uniform asymptotic expansion, we obtain expressions for 
the Bessel and Hankel functions which can readily be 
evaluated for almost arbitrary values of v and z. How- 
ever, formidable numerical cancellations still prevent us 
from using the uniform asymptotic expansions directly 
when V K, z. In this case, the region near v = z 
is bridged by the application of a convergence acceler- 
ator (Weniger transformation), applied in this case to 
the infinite asymptotic series defining the expansion for 
large v in the uniform asymptotic formulas (HU, (|18p . 
(|T^ and (pHl). The transformation is applied after the 
order v is shifted away [see Eq. (|M1) ] by a finite amount 
from the argument z using the recurrence formula satis- 
fied by the Bessel function pOl . Numerical examples are 
given in Section fill Dl 

In a typical application within physics, one needs 
Bessel functions for both real as well as complex argu- 
ment. For example, in bound-state quantum electrody- 
namics (QED), one needs to describe complex photon 
energies. Therefore, it is highly desirable to have an al- 
gorithm that is applicable in the entire complex plane, 
for a Bessel function of large order. Here, we use a com- 
bination of ideas to overcome the severe difficulties faced 
by this endeavor. We use a nonlinear sequence transfor- 
mation for the summation of asymptotic expansions for 
the Airy functions in order be able to use this expan- 
sion for small and moderate argument where the usual 
paradigm of truncating the asymptotics at the smallest 
term would otherwise yield dissatisfactory accuracy. The 
Stokes phenomenon must still be taken into account in 
terms of the complex phase of the argument of the Airy 
function. Finally, we overcome the remaining numerical 
difficulties in the summation of the uniform asymptotic 
expansion, associated with the turning point = z, by 
a simple shift of the order versus the argument, and by 
the further use of the Weniger transformation for the 
summation of the uniform asymptotic expansion in the 
transition region near v = z. 

Examples where Bessel functions of large argument 
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and order are needed, include scattering problems and 
calculations with photon propagators and fermion prop- 
agators in atomic physics and field theory, and diverse 
technical application areas. Although it is possible to re- 
cursively evaluate arrays of Bessel functions of the same 
argument with varying order, the individual evaluation of 
Bessel functions in "extreme" argument ranges remains 
an elusive problem. Here, we attempt to address this 
problem via a combination of summation algorithms, 
convergence accelerators, recurrence relations and the 
"principle of asymptotic overlap," adapted to the prob- 
lem at hand. In passing, we also address the numerically 
accurate calculation of Airy functions in the entire com- 
plex plane. We leave it to the interested reader to develop 
their own implementation of the methods described here, 
adapted to the arithmetic accuracy requirements for the 
particular application in question. 



Finally, let us indicate a few open problems in the area. 
We have outlined the general algorithm used in our eval- 
uation of Bessel functions. The description of an im- 
plementation, including an example code, possibly with 
more refined termination criteria than those indicated in 
Section IIIIBI for the Airy functions, would certainly be 
of value to the scientific community. Secondly, in view 
of the considerations outlined below in Appendix it 
seems feasible to construct an equally powerful algorithm 
based on an adaptation of the saddle point integration. 
In this case, considerable care needs to be vested into 
the turning point case « z as well. Finally, prelim- 
inary investigations (not described here in any further 
detail) indicate that the Debye expansion, which is dif- 
ferent from the uniform asymptotic expansion, may be 
used to good effect for the case of real and real ar- 
gument z, of the Bessel function. The Debye expansion 
is given in Eqs. (9.3.7),. (9.3.8), (9.3.11) and (9.3.12) of 
Ref. [12[. The Debye expansion may be combined with 
convergence accelerators. It would be instructive and 
worthwhile to construct a complementary algorithm for 
Bessel and Hankel functions valid only for real argument. 
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Appendix A: SADDLE POINT INTEGRATIONS 



Orientation 



In this Appendix, we explore alternative numerical 
procedures, with a special emphasis on numerical inte- 
gration around saddle points. Before we come to a dis- 
cussion of the procedure involved, let us briefly review 
other, alternative approaches to the calculation of Bessel 
and Hankel functions which have been discussed in the 
literature. In Ref. [l^] , it has been advocated to directly 
integrate the defining differential equation in suitable di- 
rections in the complex plane. In Ref. [26j . an analogous 
approach has been advocated for the calculation of Airy 
functions. There is a further obvious problem at the turn- 
ing point, where the argument z and order u of the Bessel 
function are nearly equal. In this case, there are huge nu- 
merical cancellations in the calculation of the expansion 
coefficients that multiply the Airy functions which are 
required in order to evaluate the uniform asymptotic ex- 
pansions. In Ref. [1^, it has been suggested to avoid 
using the expansion in inverse powers of the order of the 
Bessel function in this case, and to use a convergent ex- 
pansions in terms of a scaled complex argument Q of the 
special function. The procedure advocated in Ref. [l^ is 
applicable to the transition region v z where C ~ Oi but 
cannot be universally applied when C, is manifestly differ- 
ent from zero. Another method [29| is to expand the co- 
efficients of the asymptotic expansion into hyperasymp- 
totics. Yet, the asymptotic expansions of the coefficients 
of the asymptotic expansions (sic!) arc not applicable to 
all cases of interest and introduce another level of com- 
plexity into the problem. Hadamard series expansions as 
an alternative to the uniform asymptotics (still equal in 
limiting cases) have been investigated in Refs. |30l432l | . 
Finally, let us also mention that some special parameter 
cases of interest have received attention in the literature. 
For example, asymptotic expansions for Bessel functions 
of the third kind of imaginary order have been discussed 
in Refs. Hill. 

However, the preeminent issue in the formulation of 
possible alternative evaluation methods for Bessel and 
Hankel functions seems to be connected with a numeri- 
cal integration about saddle points in the complex plane. 
Expansions about the saddle points give rise to the above 
mentioned uniform asymptotic expansions (|12p . (|18p . 
([T^ and ([20)1 . In general, the evaluation of Hankel, and 
Bessel functions can be written in terms of two saddle 
points which dominate the paths of integration. Expand- 
ing to second order about the saddle point, we obtain 
linear contours that approximate the paths of steepest 
descent and may be used in order to evaluate the func- 
tions, approximately. This approach, which has been ad- 
vocated previously in [FH and recently in [sj , is appeal- 
ing because of its relative simplicity. Since the integrand 
of the integral defining J^{z) is in general highly oscil- 
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lating, the integration paths must be chosen with care. 
Here, we describe a somewhat simphfied approach, in 
which one approximates the ideal path by a short hnear 
segment crossing the saddle point. Our treatment is sim- 
ilar to that in |34| . however, we pay special attention to 
the case when argument and order of the Bessel function 
are almost equal, which was not considered in (33 |. 

In the case of equal argument and order, v = z, the 
phenomenon of confluence of the two saddle points needs 
to be analyzed. For v ^ the two saddle points ap- 
proach each other. For precise equality v ^ z, the two 
saddle points coalesce, and there are three paths of steep- 
est descent, complemented by three paths of steepest as- 
cent, about the saddle point. This leads to a more com- 
plicated situation with kinks in the contours of steepest 
descent, which must be used in order to evaluate the 
special functions. The investigations reported below are 
therefore of more general interest with respect to the nu- 
merical problems related to the confluence of the saddle 
points. 



2. Complex integral representations 

The defining integrals for the Hankcl functions H'iF' {z) 
and (z) read (Rez > 0, and arbitrary complex v) 

1 

7ri 



Hl^\^) = - / exp(5(i))dt, (Ala) 

-I — OO + iTT 

Hi'\z)^- exp(5(t))dt, (Alb) 

TTl Joe 



with 



g{t) = -z sinh(t) + ut = -u {y sinh(t) - t) . (A2) 

Here we have introduced the variable y — z/iy, which 
will prove convenient in the ensuing discussion. Com- 
bining the two Hankel functions to the Bessel function 
J^(z) = [Hi^\z) + Hi^\z)]/2, we obtain from Eq. (lAlall 
the integral representation for Jy{z): 



Ju{z) 



2ti\ 



— 00-]-17T 



Gxp(5(t)) dt. 



(A3) 



— oo — ITT 



The contour of integration (jASj) is arbitrary. It can be 
deformed, nevertheless, to pass through the saddle point. 
Finally, the Bessel function of the second kind is obtained 

asr,(z) = [Hi^\z)-Hl^\z)]H2\) = ^\[H^^\z)-J,{z)]. 



The definition (|Ala[) is valid provided Rez > 0, to 
ensure convergence of the integrals. The identity sinh(i — 
m) = sinh(< -I- Itt) ~ — sinh(t) is instrumental in deriving 
this condition. For Rc(z) < 0, we apply the conversion 



formulas (j6]) . The case of purely imaginary z = .t + i y 
requires another integral representation. We may use 



Jui^\y\) = cxp ( -i^TTi ) h{\y\), 



(A4) 



JA-M) = exp ( --J^TTi ) h{-\y\), (A5) 



y\ 



(A6) 
(A7) 



which expresses the Bessel functions in terms of modified 
Bessel functions Iy{x) and K^{x). For Re(a:;) > 0, these 
have integral representations 



1 r 

Iv{x) ~ — I eyip{xcos9) cos{v6)di6 
7^ Jo 



sin(:^7r) 



exp (—a; cosh t — vt) dt, 



(A8a) 



Ki,{x) = I exp (—X cosh t) cosh(z^t)dt . (A8b) 

"'0 

Negative arguments can be transformed to positive argu- 
ments via the relations 

Iu{—x) = exp(z^7ri) /^(x), (A9a) 
K^{-x) = exp(— I'Tri) Ky{x) — TT\Iy{x). (A9b) 



3. Paths of steepest descent 



We recall the definition of the integrand in the integral 
representation of the Bessel function. 



g{t) = -V {y sinh(i) - t) . 



(AlO) 



Expanded to second order about the saddle point t ^ tj, 
to be specified below, one obtains 

9{t) = 9{tj) + 9'{tj) {t - tj) + ig"(i,) {t - tjf + ... 
= — V (ysinh(tj) — tj) + v {\ — y cosh(tj)) (t — tj) 



-vy cosh{t,){t-tjf + 



(All) 



The method of numerical evaluation of the functions 
Ju{z) and Yi,{z) amounts to finding approximations to 
the integrals (jAlaj) . ()Alb[) by quadrature. The starting 
point is to find the saddle points t± of the integrand in 
[3|) . They satisfy g'{t±) = = ycosht± — 1, and are 
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given as 



t+ = In 



y 



= a+ + i 6+, 



(Al2a) 



t_ = \u\^ ^ j ^ + i b- = , (A12b) 

for |Imt±| < TT, and the notation t± ^ a± + ib± (real 
o-ij b±) is introduced for later use. Multiplying the ar- 
guments of the logarithms in Eq. (jA12p . we sec that the 
saddle points fulfill the relation f+ + t_ = 0. Further- 
more, we have Ret_ < (and Rei+ > 0). For the case 
of equal argument and order, the two saddle point coa- 
lesce at the origin, 



i+ = t_ = 0. 



(A13) 



The path of steepest descent (PSD) is a path t in the 
complex plane, passing through the saddle point, such 
that the imaginary part of the exponent git) is constant 
along this path. This means that the integrand exp {g(t)) 
does not oscillate. Of course, there are further saddle 
points in the complex plane, described by the formula 
t± + 2n7ri, n integer, but these are not needed in our 
investigation. 

In addition, the amplitude of the integrand exp {g{t)) 
decreases monotonically as we go along the path away 
from the saddle point. Together, these properties make 
numerical quadrature along the PSD easy. In order find 
an expression for the PSD, we write 



(A14) 



and require that the imaginary part of g{t) remains con- 
stant along the path of integration. We write z = zr+I zj, 
with zr = Rez and zi = Imz, and in complete analogy, 
V = i/n + ivi. Then, 

lmg{t±) =lmg{t) (A15) 
= ~zi sinh ^ cos 77 — cosh ^ sin i] + vi + rj. 



In general, Eq. (|A15[) cannot be solved analytically. The 
limiting behavior as |^| — 00 may however be deduced. 
Namely, for |^| — >• 00, one has |sinh^| |^|, |cosh^| ^ 
1^1, and tanh^ — >• 1 for ^ — > ±00. Therefore, the asymp- 
totic solutions are given by 

tan?7 = — = tan(argz) , ^ —00, (Al6a) 

tan 77= ^- = — tan (arg z) , ^ — > 00. (A16b) 



A possible approach now is to solve Eq. (|A15p numeri- 
cally, and subsequently integrate along the numerically 
obtained PSD [SJ] . A helpful discussion of numerical as- 
pects related to contour integrals in the complex plane is 
given in Chapter 5 Ref. [l8| . 



As we shall see, to obtain modest accuracy it is enough 
to integrate along a short linear segment F, which ap- 
proximates the true PSD close to the saddle point. We 
have previously defined the saddle points tj {j = ±) and 
their real and imaginary parts a j and bj . Linear approx- 
imations to the paths of steepest descent in the complex 
t plane are obtained as follows. 



^ = i-j iv) = iv- bj) Kj +aj +i 77, = ± , (A17a) 
tj = tj {bj) = % + i bj , ^l^Ml ^Kj+i, (A17b) 



where Kj is real. That is, we parameterize the integral 
in (jA3|) (similarly for the Yy{z) function) as a function of 
rj = Im t as 



hi. 27ri Jr, 



— / exp(g(tj(r/))) d?7, 



E 



di,(77) 



(A18) 

where the Tj arc line segments along the paths that define 
linear approximations to the contours of steepest descent, 
with j e {-f-, — }. The values of the parameters rymin and 
f/max are chosen so that a specific final accuracy is reached 
for the approximation to the integral. 



The sum in (|A18p may run over one, or both of the 
saddle points t± , depending on the saddle point configu- 
ration. Our expansion to second order about the saddle 
point ensures that no oscillations occur in the integrand 
up to this order. Therefore, no imaginary part will be 
incurred in the integrand upon expansion about the sad- 
dle point to second order in t — tj. However, if one goes 
beyond second order, oscillations will occur upon using 
the linear approximation (jA17p . 



The coefficients Kj in (|A17p can be found by expand- 
ing Eq. (jA15|) to second order around the saddle point 
tj , which yields 

= lm [z(t - tj)'^ fim\i{tj)\ 
= Im {z [(C - a,)2 _ (,y - bj)^ + 2i(e - a,){7^ - b,)] 
X [sinh aj cos bj + i cosh aj sin foj] } . (A19) 



Solving for ^ as a function of ?/. we obtain 



{V~b,) ± ^Bj + 1) 



(A20) 



where 



B, 



zr sinh a j cos bj — zi cosh aj sin bj 
zr cosh aj sin bj + zj sinh aj cos bj 
cot arg(z sinh tj). (A-21) 
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It will be discussed below how to estimate the upper and 
lower integration limits rj^-m and ?7max- The coefficient 
Kj is given by 

K, = -B,±^fBf7l, (A22) 

and the ± sign needs to be fixed. Indeed, there are two 
solutions for each saddle point, one of which corresponds 
to the PSD. The other solution yields the path of steepest 
ascent, along which the integrand exp{g{t)) increases. 
The condition 

* arg(zsinhij) G (-7r,7r), (A23a) 

{Kj + if = ^K] + 1 exp (-i , (A23b) 

ensures that the term of order [t — tjf cx rf in Eq. (|Alip 
leads to an exponential decrease (instead of increase) of 
the integrand around the saddle point. It is fulfilled if we 
choose 

K, = -B,- ^B] + 1 (A24a) 
= - -y/cot2 (i^-) + lexp(-i*) , «'>0, 

= -Bj + ^B] + 1 (A24b) 
= \/cot2(i5') +lexp(-i«') , *<0, 

depending on the sign of 4*. 

The case v = z needs to be considered separately. The 
saddle point coalesce at i+ =t- =0. The second-order 
terms vanishes (sinhtj = 0), and the linear approxima- 
tion to the PSD is obtained by expanding (|A15[) to third 
order. In this case, the saddle point becomes a "triple 
saddle point" which instead of a cross formation has the 
shape of a double saddle point, or of a "six- fold star" 
with three paths of steepest descent, and three paths of 
steepest ascent in between. The corresponding term from 
Eq. (jAlip in this case is 

- ^ cosh(i, = 0) (i - t,f = (A25) 

because, for the confluence of the two saddle points, we 
have V = z and tj = 0. The saddle point equation thus 
simplifies and reads 

Im(zt3) = 0. (A26) 
This has three solutions, 

tk = 77(cot ifk + i) = -. exp(i(^fe), (A27) 

sm ipk 

where ipk = {kn — argz)/3, and k e {0,1,2}, and the 
last equality indicates that the confluent PSDs can be pa- 
rameterized conveniently in terms oi s = rj/ sin ip^ . The 
PSDs derived in Eq. (|A27|) can be used not only when 



1/ = z, but also in the case where v and z are close, but 
not exactly equal. The three solutions given in Eq. (|A27p 
change from paths of steepest ascent to paths of steepest 
descent, depending on the sign of 77. In particular, when 
argz = 0, then the PSD for Jy(z) close to the saddle 
point t ^ satisfies argt = —2tt/3 when Imt < and 
argt = 27r/3 when Imt > 0. The PSD thus approaches 
the origin from the lower left and departs again to the 
upper left [see also Fig. |6ljc)]. 

The confluence of the two saddle points in the case v = 
z represents the major obstacle in a numerical treatment 
of the Bessel function, in the context of the saddle-point 
integration. 



4. Numerical integration 

In terms of the PSD configurations, there are three 
cases that should be considered separately, depending on 
the sign of Im (<:+), where is the saddle point with 
Re(t+) > 0. The saddle point configuration depends on 
this sign (33 |. 

Case 1: If Invit^) < and consequently Im(t_) > 0, 
the configuration of the saddle points and PSDs in this 
case is illustrated in Fig. [6lja). As is gathered from this 
plot, one saddle point contributes to the value of J^{z), 

and the other to the value of H^\z). The remain- 
ing function Yi,{z) can then be computed as Y^{z) ~ 

-i[Hi^\z) - Juiz)]. Case 2: If Im(t_) < and there- 
fore Im (t+) > 0, the situation is as shown in Fig. [6l^b). 
To obtain Ji/(z), one has to pass both saddle points, and 
thus pick up a contribution from both. Case 3: The case 
Imt+ = Imt_ = actually needs a separate considera- 
tion. It is shown in Fig. |6l[c), but only in the special case 
of coalescing saddle points (see the discussion below) . In 
general, the saddle points may be separated, still, with 
Imt_|_ = Imt_ = 0. 

For V K, z^ but not exact equality, the two saddle points 
almost coalesce. In Fig. HKc), we investigate the case 

V — z^ with = t_ = 0. In this case, the saddle point 
has merged together into a double saddle point at the 
origin, and Jv{z) and Hi ' {z) are calculated by following 
the appropriate paths indicated by red and blue arrows, 
respectively. The linear segment approximation to the 
true PSD [dashed lines in Fig. HJc)] are seen to be valid 
rather far from the saddle point, which implies that the 
linear approximation to the PSD may be used to good 
effect. The given linear approximation to the PSD for 

V = z may be used also for the case 1 1 — 7 1 < h with 
given h. In practice, the value of h increases with increas- 
ing |z|, for a given prescribed accuracy of the numerical 
evaluation. As a rule of thumb, the imaginary part of 
git) along the effective path should be of order 1, where 
the length of the path is estimated by Eq. (|A34[) . 
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FIG. 6: (Color online.) PSD for different values of the ar- 
gument z and order v. In (a), f = (1000 + l/5)e'''/^°, 
z = (1200 + 6/25)e2'"/^ in (b), u = (1000 + l/5)e'"/^ 
2 = (1100+ll/50)e'''/^", and in (c), z^u = (1000+l/5)e'''''^ 
Red lines show PSDs passing through the saddle points (black 
dots), blue lines paths of steepest ascent passing through the 
saddle points. For orientation purposes, thin black lines are 
drawn at Imt — ±7r. Red arrows indicate the integration 
path to obtain Ji,{z), and blue arrows the paths for H^^\z) 
[(a) and (b)] and H^'^\z) [(c)]. In (a) and (b), the green curves 
are curves satisfying Eq. (|A15|I . but not passing through the 
saddle point. The slopes of the linear segments (|A20|) . (IA27|) 
used in the numerical evaluation are shown with dashed lines. 
To obtain 10 significant figures, one needs to perform the inte- 
gration along a linear segment of length approximately equal 
to the size of the black dots. 



To gain a little more insight into the difficulties asso- 
ciated with the saddle points approaching each other, we 
discuss the transition from Fig.[6l[a) to|6l[c), by a suitable 
change of the parameters. The PSD for Jv{z) (red ar- 
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FIG. 7; (Color online.) Real (dashed, blue line) and imag- 
inary (red, solid line) parts of the integrand F(r\) from 
Eq. The parameter values are v = (500 -I- l/5)e^'''/^'', 

and z = (500 + l/10)e''''^/^'^+^°"'''. The limits of the 77-axis 
correspond to those obtained from the estimate (|A34|l 



rows) then changes from having smooth curvature at the 
saddle point [Fig. El^a)] to having a "kink" at the saddle 
point [Fig.[6l^c)]. Therefore, even if the two saddle points 
are separated, at some point the linear approximation to 
the PSD will only hold for a very short segment of the 
true PSD, and is unsuitable. For some small distance 
between the saddles, it is instead better to pretend that 
the suitable contour is given by the double saddle point 
and to use the dashed lines in Fig. IH^c), even if the saddle 
points do not exactly coalesce. We illustrate the behav- 
ior of the integrand along the approximate PSD when 
v fa z in Fig. [T] Shown are real and imaginary parts of 
the integrand F{ri), where 



J. (2) 



Fi,) = -e. 



aitin)) 



drj 



(A28) 



(A29) 



and [see Eq. (|A27P ] 



cot 



TT — arg z 



27r — arg z 
cot I + 1 



if 77 < , 
, if 77 > 0. 



(A30) 



In Fig. [71 we have h 



1 - - 



3 X 10 



-3 



The integration limits for the numerical calculations 
can be chosen as follows. The primary objective is to 
estimate how far out from the saddle point the integra- 
tion should be cut off, that is to estimate the limits rj^nn 
and ?/max in Eq. (|A18I) . In the case of separate sad- 
dle points, we expand g{t) around the saddle point up to 
second order, to obtain an integral of the form 



hix) 



-Cis' 



ds , 



(A31) 



17 



where X and Ci are positive real constants. It follows 
that if we want relative precision e = |1 — /i(X)//i(oo)| 
we must cut the integral at 



X 



erfc'^e) 



(A32) 



where crfc~^(-) is the inverse of the complementary error 
function. For coalescing saddle points, also the second 
derivative of g(t) vanishes, and we must instead consider 
integrals like 



I2{X)=^ f e 



ds. 



Here, we should cut at 
X = 



C2 



(A33) 



(A34) 



to obtain a prescribed precision e. In Eq. (|A34[) . r^^(-, •) 
is the inverse of the incomplete gamma function with re- 
spect to its second argument. These estimates assumes 
a monotonically decreasing integrand along the path of 
steepest descent, and it is clear that it is impossible 
to reach arbitrary precision due to incurred oscillations 
along the linear line segments approximating the paths 
of steepest descent, as one travels too far from the saddle 
point. 

In Ref. [13], the authors also advocate to use an ap- 
proximation of the saddle point contour by linear line 
segments. They first observe that it is possible to solve 



the PSD equation (|A15|) numerically, but then show that 
if oscillations do not induce prohibitive numerical insta- 
bility, equally accurate results may be obtained by just 
taking line segments. Because it is time-consuming to 
calculate the PSD numerically, they conclude that the 
line-segment approach should be favored. However, as 
mentioned previously, this approach necessarily breaks 
down at some level of precision, especially for large val- 
ues of |z|. To reach (in principle) arbitrary precision for 
large magnitudes of the order and argument with the 
saddle point method, it seems unavoidable to follow the 
true PSD as closely as possible. 



An optimized, hybrid algorithm should perform the 
PSD stepping and the integration in parallel: First, one 
advances on the PSD by one step, computes the inte- 
gral on the segment connecting the previous point by 
the new point, checks if the accuracy demand is sat- 
isfied, if not, one takes another step. The step length 
could be adjusted according to how close the condition 
lm{g{t±)) = lm{g{t)) is satisfied. Since the integrand 
decreases exponentially along the PSD, and does not os- 
cillate, it should be straightforward to estimate when to 
stop the stepping/integration. Still, for this algorithm 
to be universally applicable, it might be necessary to in- 
vestigate the overlap regions where the saddle-point con- 
figurations versus the double saddle point configurations 
should be used, and to develop special routines that deal 
with the line segments joining the two saddle points in 
overlapping regions. We have not pursued this endeavor 
any further in the current work and leave it as an open 
problem. 
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